%% Output Figures for Baseline Solution
for( option_Q = [1]    )
MainPath = '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model/Codes_2017-09-15';

if( option_Q==1 )
    cd( [MainPath, '/Results-Constant Q'] );
elseif( option_Q==2 )
    cd( [MainPath, '/Results-Linear Q'] );
elseif( option_Q ==3 )
    cd( [MainPath, '/Results-Quadratic Q'] );
else
    cd( [MainPath, '/Results-Pivot'] );
end

load( [ 'solution(par set 1).mat' ] )
% Recovery a sponge Illustration.
h = figure('position', [0, 0, 400, 300]);
sub_w_index  =  w_grid>0.01 & w_grid<0.3 ;   subindex_ode = w_ode45>0.01 & w_ode45<0.3 ;
i1 = 14;  i2 = 8; 
plot( w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i1) ,  w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i2),  '--'  , w_ode45(subindex_ode), muw_ratio_ode45(subindex_ode), ':m',   w_ode45(subindex_ode),  0*w_ode45(subindex_ode), '-.k',   'linewidth', linewidth );    
title('growth rate of w'); xlabel('w');      legend(   ['Bank-run Model with \lambda=', num2str(lambda_grid(i1))] ,   ['Bank-run Model with \lambda=', num2str(lambda_grid(i2))],  'Model without Bank-run'  );
saveTightFigure( h,  [ 'Recovery_As_Sponge.pdf' ]  );  close all;

load( [ 'solution(par set 2).mat' ] )
h = figure('position', [0, 0, 900, 300]);
linewidth=1.5;  reset(h);  
i1 = 5;      i2 = 15 ;  i3 = 25; subindex_lambda = lambda_grid<0.5;
subplot(1,3,1);  plot( [0;lambda_grid],  [xK_ode45(i1), xK_matrix(i1,:)]  ,  [0;lambda_grid], [xK_ode45(i2),xK_matrix(i2,:)],  '--'  ,[0;lambda_grid] ,  [xK_ode45(i3),xK_matrix(i3,:)]  , ':',    'linewidth', linewidth );    title('banker capital holding x^K'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(1,3,2);  plot( [0;lambda_grid(subindex_lambda) ],   [ Rp_ode45(i1), RpH_matrix(i1,subindex_lambda) ] ,  [0;lambda_grid(subindex_lambda)], [ Rp_ode45(i2),RpH_matrix(i2,subindex_lambda)],  '--'  ,[0;lambda_grid(subindex_lambda)],   [ Rp_ode45(i3),RpH_matrix(i3,subindex_lambda)] , ':',   'linewidth', linewidth);    title('risk premium'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(1,3,3);  plot( [0;lambda_grid],  [ 0, kappap_matrix(i1,:) ] ,  [0;lambda_grid], [ 0,kappap_matrix(i2,:)],   '--'  ,[0;lambda_grid] ,  [0,kappap_matrix(i3,:)]  , ':' ,   'linewidth', linewidth );    title('price jump \kappa^p'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
saveTightFigure( h, 'Pre-Crisis.pdf' );
close all;


%% Output Figures for Main Solutions
for( parset_choice = 1:3 )

clearvars -except parset_choice  MainPath  % Note: it is important to clear all other variables, and only leave parset_choice.
load( [ 'solution(par set ', num2str(parset_choice)  ,').mat' ] )

% Output error trend
h = figure; plot( log(errors) );  xlabel('iteration'); ylabel('log(absolute error)')
saveTightFigure( h,  ['ErrorTrend(par set ', num2str(parset_choice)  ,').pdf']  );  close all;

% Output Figure 1
clf('reset');
h = figure('position', [0, 0, 900, 920]);
linewidth=1.5;  reset(h);  
i1 = 15;  i2 = 7; i3=2;
subplot(4,3,1);  plot(  w_grid, p_2D_new(:,i1),   w_grid, p_2D_new(:,i2),  '--'  ,  w_grid, p_2D_new(:,i3),  '-.'  , w_ode45 , p_ode45 , ':m',   'linewidth', linewidth  );    title('price p');     xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], ['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
sub_w_index  =  w_grid>0.02 ;   subindex_ode = w_ode45>0.02;
subplot(4,3,2);  plot(  w_grid(sub_w_index), pw_matrix(sub_w_index, i1),  w_grid(sub_w_index), pw_matrix(sub_w_index, i2),  '--'  ,   w_grid(sub_w_index), pw_matrix(sub_w_index, i3),  '-.' , w_ode45(subindex_ode) , pw_ode45(subindex_ode) , ':m',   'linewidth', linewidth  );    title('price derivative p_w');  xlabel('w');   legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], ['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,3);  plot(  w_grid(sub_w_index), sigmap_matrix(sub_w_index, i1),  w_grid(sub_w_index), sigmap_matrix(sub_w_index, i2),  '--'  ,  w_grid(sub_w_index), sigmap_matrix(sub_w_index, i3),  '-.' , w_ode45(subindex_ode) , sigmap_ode45(subindex_ode) , ':m',  'linewidth', linewidth );    title('price volatility \sigma^p');  xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], ['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subindex_ode = xK_ode45<15& (psi_ode45<=i1) & w_ode45~=0 & w_ode45~=1;
subplot(4,3,4);  plot( w_grid(sub_w_index), kappap_matrix(sub_w_index, i1),  w_grid(sub_w_index), kappap_matrix(sub_w_index, i2),   '--'  , w_grid(sub_w_index), kappap_matrix(sub_w_index, i3),   '-.' ,  w_ode45 ,  0*N_ode45  , ':m' ,   'linewidth', linewidth );    title('price jump \kappa^p'); xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))],  '\lambda=0'  );
subplot(4,3,5);  plot( w_grid(sub_w_index), xK_matrix(sub_w_index, i1),  w_grid(sub_w_index), xK_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), xK_matrix(sub_w_index, i3),  '-.'  ,  w_ode45(subindex_ode) ,  xK_ode45(subindex_ode)  , ':m',    'linewidth', linewidth );    title('banker capital holding x^K'); xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,6);  plot( w_grid(sub_w_index), xg_matrix(sub_w_index, i1),  w_grid(sub_w_index), xg_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), xg_matrix(sub_w_index, i3),  '-.'  ,   w_ode45(subindex_ode) ,  xg_ode45(subindex_ode)  ,  ':m',    'linewidth', linewidth);    title('banker liquidity holding x^g'); xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,7);  plot( w_grid(sub_w_index), liq_prem_matrix(sub_w_index, i1),  w_grid(sub_w_index), liq_prem_matrix(sub_w_index, i2),  '--'  ,  w_grid(sub_w_index), liq_prem_matrix(sub_w_index, i3),  '-.' ,  w_ode45,  zeros(N_ode45,i1), ':m',  'linewidth', linewidth);    title('liquidity premium'); xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], ['\lambda=', num2str(lambda_grid(i3))],'\lambda=0'  );
subplot(4,3,8);  plot( w_grid(sub_w_index), RpH_matrix(sub_w_index, i1),  w_grid(sub_w_index), RpH_matrix(sub_w_index, i2),  '--', w_grid(sub_w_index), RpH_matrix(sub_w_index, i3),  '-.',    w_ode45,   muR_ode45 + AH./p_ode45 , ':m',   'linewidth', linewidth);    title('risk premium'); xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,9);  plot( w_grid(sub_w_index), SharpeRatioH_matrix(sub_w_index, i1),  w_grid(sub_w_index), SharpeRatioH_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), SharpeRatioH_matrix(sub_w_index, i3),  '-.' ,  w_ode45(sub_w_index), SharpeRatio_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('Sharpe ratio'); xlabel('w');   legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,10);  plot( w_grid(sub_w_index), SR_adjusted_matrix(sub_w_index, i1),  w_grid(sub_w_index), SR_adjusted_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), SR_adjusted_matrix(sub_w_index, i3),  '-.' ,  w_ode45(sub_w_index), SharpeRatio_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('adjusted Sharpe ratio'); xlabel('w');   legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subplot(4,3,11);  plot( w_grid(sub_w_index), r_matrix(sub_w_index, i1),  w_grid(sub_w_index), r_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), r_matrix(sub_w_index, i3),  '-.' ,  w_ode45(sub_w_index), r_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('interest rate'); xlabel('w');   legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
subindex_ode  =  muw_ode45<0.3 ;
subplot(4,3,12);  plot( w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i1),  w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i2),  '--'  , w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i3),  '-.' ,  w_ode45(subindex_ode), muw_ratio_ode45(subindex_ode),  ':m',  'linewidth', linewidth);    title('growth rate of w'); xlabel('w');   legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],['\lambda=', num2str(lambda_grid(i3))], '\lambda=0'  );
saveTightFigure( h,  [ 'Baseline_w(par set ', num2str(parset_choice)  ,').pdf' ]  );  close all;


% Output Figure 2
clf('reset');
h = figure('position', [0, 0, 900, 920]);
linewidth=1.5;  reset(h);  
i1 = 5;      i2 = 15 ;  i3 = 25;
subplot(4,3,1);  plot(  [0;lambda_grid] , [ p_ode45(i1), p_2D_new(i1,:) ],   [0;lambda_grid],  [ p_ode45(i2), p_2D_new(i2,:)],  '--'  ,  [0;lambda_grid] , [ p_ode45(i3), p_2D_new(i3,:)] , ':',   'linewidth', linewidth  );    title('price p');     xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))],  ['w=', num2str(w_grid(i3))]  );
subplot(4,3,2);  plot(  [lambda_grid],  [ pw_matrix(i1,:)] ,  [lambda_grid],  [ pw_matrix(i2,:)],  '--'  ,[lambda_grid] ,  [pw_matrix(i3,:)] , ':',   'linewidth', linewidth  );    title('price derivative p_w');  xlabel('lambda');   legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
%subplot(4,3,3);  plot(  [0;lambda_grid],  [sigmap_ode45(i1),sigmap_matrix(i1,:)]  ,    [0;lambda_grid],  [sigmap_ode45(i2),sigmap_matrix(i2,:)],  '--'  , [0;lambda_grid],  [sigmap_ode45(i3),sigmap_matrix(i3,:)], ':',  'linewidth', linewidth );    title('price volatility \sigma^p');  xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,3);  plot(  [lambda_grid],  [sigmap_matrix(i1,:)]  ,    [lambda_grid],  [sigmap_matrix(i2,:)],  '--'  , [lambda_grid],  [sigmap_matrix(i3,:)], ':',  'linewidth', linewidth );    title('price volatility \sigma^p');  xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,4);  plot( [0;lambda_grid],  [ 0, kappap_matrix(i1,:) ] ,  [0;lambda_grid], [ 0,kappap_matrix(i2,:)],   '--'  ,[0;lambda_grid] ,  [0,kappap_matrix(i3,:)]  , ':' ,   'linewidth', linewidth );    title('price jump \kappa^p'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,5);  plot( [0;lambda_grid],  [xK_ode45(i1), xK_matrix(i1,:)]  ,  [0;lambda_grid], [xK_ode45(i2),xK_matrix(i2,:)],  '--'  ,[0;lambda_grid] ,  [xK_ode45(i3),xK_matrix(i3,:)]  , ':',    'linewidth', linewidth );    title('banker capital holding x^K'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,6);  plot( [0;lambda_grid],  [xg_ode45(i1), xg_matrix(i1,:)]  ,  [0;lambda_grid], [xg_ode45(i2),xg_matrix(i2,:)],  '--'  ,[0;lambda_grid] ,   [xg_ode45(i3),xg_matrix(i3,:) ] ,  ':',    'linewidth', linewidth);    title('banker liquidity holding x^g'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,7);  plot( [0;lambda_grid],   [0, liq_prem_matrix(i1,:)] ,  [0;lambda_grid],  [0, liq_prem_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [0,  liq_prem_matrix(i3,:)], ':',  'linewidth', linewidth);    title('liquidity premium'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,8);  plot( [0;lambda_grid],   [ Rp_ode45(i1), RpH_matrix(i1,:) ] ,  [0;lambda_grid], [ Rp_ode45(i2),RpH_matrix(i2,:)],  '--'  ,[0;lambda_grid],   [ Rp_ode45(i3),RpH_matrix(i3,:)] , ':',   'linewidth', linewidth);    title('risk premium'); xlabel('lambda');  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,9);  plot( [0;lambda_grid], [SharpeRatio_ode45(i1), SharpeRatioH_matrix(i1,:)] ,  [0;lambda_grid], [SharpeRatio_ode45(i2),SharpeRatioH_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [SharpeRatio_ode45(i3),SharpeRatioH_matrix(i3,:)] ,  'linewidth', linewidth);    title('Sharpe ratio'); xlabel('lambda');   legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,10);  plot( [0;lambda_grid], [SharpeRatio_ode45(i1), SR_adjusted_matrix(i1,:)] ,  [0;lambda_grid], [SharpeRatio_ode45(i2),SR_adjusted_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [SharpeRatio_ode45(i3),SR_adjusted_matrix(i3,:)] ,':', 'linewidth', linewidth);    title('adjusted Sharpe ratio'); xlabel('lambda');   legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,11);  plot( [0;lambda_grid], [r_ode45(i1), r_matrix(i1,:)] ,  [0;lambda_grid], [r_ode45(i2),r_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [r_ode45(i3),r_matrix(i3,:)] , ':', 'linewidth', linewidth);    title('interest rate'); xlabel('lambda');   legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(4,3,12);  plot( [ lambda_grid], [  muw_ratio_matrix(i1,:)] ,  [ lambda_grid],  [ muw_ratio_matrix(i2,:)],  '--'  ,[ lambda_grid],  [muw_ratio_matrix(i3,:)] ,':',  'linewidth', linewidth);    title('growth rate of w'); xlabel('lambda');   legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
saveTightFigure( h,  [ 'Baseline_lambda(par set ', num2str(parset_choice)  ,').pdf' ]  );  close all;

end


%% Compare Different Q
load('solution(par set 1).mat'); 
S1=load('solution(par set 1).mat');     S2=load('solution(par set 3).mat');  
labelS1_1 = ['Q=', num2str(S1.multiplier)];  labelS1_2= ['Q=', num2str(S1.multiplier), ' no jump'];    
labelS2_1 = ['Q=', num2str(S2.multiplier)];  labelS2_2= ['Q=', num2str(S2.multiplier), ' no jump'];    
if( option_Q == 4 )
    labelS1_1 = ['Countercyclical Q'];  labelS1_2= ['Countercyclical Q, no jump'];    
    labelS2_1 = ['Cyclical Q'];  labelS2_2= ['Cyclical Q, no jump'];  
end

%% ****** w dimension ******
for( index_lambda = [4 11]  )
index_w  =  w_grid>0.02 ;   index_ode=w_ode45>0.02;  linewidth=1;
h = figure('position', [0, 0, 900, 920]);  
% Part 1: price comparison
subplot(4,3,1);
plot(  S1.w_grid , S1.p_2D_new(:,index_lambda),  '--', 'linewidth', linewidth  );    title('price p');     xlabel('w');  hold on;
plot(  S2.w_grid , S2.p_2D_new(:,index_lambda),  'linewidth', linewidth  );    title('price p');     xlabel('w'); 
% Part 2: derivative comparison
subplot(4,3,2);
plot(  S1.w_grid(index_w) , S1.pw_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );    title('price derivative p_w');     xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.pw_matrix(index_w,index_lambda),  'linewidth', linewidth  );    title('price derivative p_w');     xlabel('w'); 
legend( labelS1_1, labelS2_1 ); 
% Part 3: sigmap comparison
subplot(4,3,3);
plot(  S1.w_grid(index_w) , S1.sigmap_matrix(index_w,index_lambda), '--', 'linewidth', linewidth  );     title('price volatility \sigma^p');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.sigmap_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('price volatility \sigma^p');     xlabel('w'); 
% Part 4: kappap comparison
subplot(4,3,4);
plot(  S1.w_grid(index_w) , S1.kappap_matrix(index_w,index_lambda), '--',   'linewidth', linewidth  );     title('price jump \kappa^p');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.kappap_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('price jump \kappa^p');     xlabel('w'); 
% Part 5: xK comparison
subplot(4,3,5);
plot(  S1.w_grid(index_w) , S1.xK_matrix(index_w,index_lambda), '--', 'linewidth', linewidth  );     title('banker capital holding x^K');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.xK_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('banker capital holding x^K');     xlabel('w'); 
% Part 6: xg comparison
subplot(4,3,6);
plot(  S1.w_grid(index_w) , S1.xg_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('banker liquidity holding x^g');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.xg_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('banker liquidity holding x^g');     xlabel('w'); 
% Part 7: liquidity premium comparison
subplot(4,3,7);
plot(  S1.w_grid(index_w) , S1.liq_prem_matrix(index_w,index_lambda) , '--',  'linewidth', linewidth  );     title('liquidity premium');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.liq_prem_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('liquidity premium');     xlabel('w'); 
% Part 8: risk premium comparison
subplot(4,3,8);
plot(  S1.w_grid(index_w) , S1.RpH_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('risk premium');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.RpH_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('risk premium');     xlabel('w'); 
% Part 9: Sharpe Ratio comparison
subplot(4,3,9);
plot(  S1.w_grid(index_w) , S1.SharpeRatioH_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('Sharpe Ratio');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.SharpeRatioH_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('Sharpe Ratio');     xlabel('w'); 
% Part 10: adjusted Sharpe Ratio comparison
subplot(4,3,10);
plot(  S1.w_grid(index_w) , S1.SR_adjusted_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('adjusted Sharpe Ratio');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.SR_adjusted_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('adjusted Sharpe Ratio');     xlabel('w'); 
% Part 11: interet rate comparison
subplot(4,3,11);
plot(  S1.w_grid(index_w) , S1.r_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('interest rate');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.r_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('interest rate');     xlabel('w'); 
% Part 12 : growth rate comparison
subplot(4,3,12);  index_ode = w_ode45>0.01;
plot(  S1.w_grid(index_w) , S1.muw_ratio_matrix(index_w,index_lambda), '--', 'linewidth', linewidth  );     title('growth rate of w');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.muw_ratio_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('growth rate of w');     xlabel('w'); 

saveTightFigure( h,  [ 'LiquidityEffects_w(lambda=',  num2str(  S1.lambda_grid( index_lambda )  )  ,').pdf' ]  );  close all;
end

%% lambda dimension
for(  w_for_lambda = [6  27  40]  )
clf('reset');  h = figure('position', [0, 0, 900, 920]);  
% Part 1: price comparison
subplot(4,3,1);
plot(  [0;S1.lambda_grid] , [ S1.p_ode45(w_for_lambda), S1.p_2D_new(w_for_lambda,:)]  ,  '--', 'linewidth', linewidth  );    title('price p');     xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [ S2.p_ode45(w_for_lambda), S2.p_2D_new(w_for_lambda,:)],  'linewidth', linewidth  );    title('price p');     xlabel('\lambda'); 
% Part 2: derivative comparison
subplot(4,3,2);
plot(  [0;S1.lambda_grid] , [S1.pw_ode45(w_for_lambda),S1.pw_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );    title('price derivative p_w');     xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.pw_ode45(w_for_lambda),S2.pw_matrix(w_for_lambda,:)], 'linewidth', linewidth  );    title('price derivative p_w');     xlabel('\lambda'); 
legend( labelS1_1, labelS2_1 ); 
% Part 3: sigmap comparison
subplot(4,3,3);
plot(  [0;S1.lambda_grid] ,  [S1.sigmap_ode45(w_for_lambda),S1.sigmap_matrix(w_for_lambda,:)] , '--', 'linewidth', linewidth  );     title('price volatility \sigma^p');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] ,  [S2.sigmap_ode45(w_for_lambda),S2.sigmap_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('price volatility \sigma^p');     xlabel('\lambda'); 
% Part 4: sigmap comparison
subplot(4,3,4);
plot(  [0;S1.lambda_grid] , [0,S1.kappap_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('price jump \kappa^p');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [0,S2.kappap_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('price jump \kappa^p');     xlabel('\lambda'); 
% Part 5: xK comparison
subplot(4,3,5);
plot(  [0;S1.lambda_grid] , [S1.xK_ode45(w_for_lambda),S1.xK_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('banker capital holding x^K');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.xK_ode45(w_for_lambda),S2.xK_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('banker capital holding x^K');     xlabel('\lambda'); 
% Part 6: xg comparison
subplot(4,3,6);
plot(  [0;S1.lambda_grid] , [S1.xg_ode45(w_for_lambda),S1.xg_matrix(w_for_lambda,:)], '--', 'linewidth', linewidth  );     title('banker liquidity holding x^g');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.xg_ode45(w_for_lambda),S2.xg_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('banker liquidity holding x^g');     xlabel('\lambda'); 
% Part 7: liquidity premium comparison
subplot(4,3,7);
plot(  [0;S1.lambda_grid] , [0,S1.liq_prem_matrix(w_for_lambda,:)] , '--',  'linewidth', linewidth  );     title('liquidity premium');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [0,S2.liq_prem_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('liquidity premium');     xlabel('\lambda'); 
% Part 8: risk premium comparison
subplot(4,3,8);
plot(  [0;S1.lambda_grid] , [S1.Rp_ode45(w_for_lambda),S1.RpH_matrix(w_for_lambda,:)], '--', 'linewidth', linewidth  );     title('risk premium');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.Rp_ode45(w_for_lambda),S2.RpH_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('risk premium');     xlabel('\lambda'); 
% Part 9: Sharpe Ratio comparison
subplot(4,3,9);
plot(  [0;S1.lambda_grid] , [S1.SharpeRatio_ode45(w_for_lambda),S1.SharpeRatioH_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('Sharpe Ratio');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.SharpeRatio_ode45(w_for_lambda),S2.SharpeRatioH_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('Sharpe Ratio');     xlabel('\lambda'); 
% Part 10: adjusted Sharpe Ratio comparison
subplot(4,3,10);
plot(  [0;S1.lambda_grid] , [S1.SharpeRatio_ode45(w_for_lambda),S1.SR_adjusted_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('adjusted Sharpe Ratio');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.SharpeRatio_ode45(w_for_lambda),S2.SR_adjusted_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('adjusted Sharpe Ratio');     xlabel('\lambda'); 
% Part 11: interest rate
subplot(4,3,11);
plot(  [0;S1.lambda_grid] , [S1.r_ode45(w_for_lambda),S1.r_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('interest rate');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.r_ode45(w_for_lambda),S2.r_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('interest rate');     xlabel('\lambda'); 
% Part 12: growth rate of w
subplot(4,3,12);
%plot(  [0;S1.lambda_grid] , [S1.muw_ratio_ode45(w_for_lambda),S1.muw_ratio_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('growth rate of w');  xlabel('\lambda');  hold on;
plot(  [S1.lambda_grid] , [S1.muw_ratio_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('growth rate of w');  xlabel('\lambda');  hold on;
%plot(  [0;S1.lambda_grid] , [S1.muw_ratio_ode45(w_for_lambda),S1.muw_ratio_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('growth rate of w');     xlabel('\lambda'); 
plot(  [S2.lambda_grid] , [S2.muw_ratio_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('growth rate of w');     xlabel('\lambda'); 

saveTightFigure( h,  [ 'LiquidityEffects_lambda(w=',  num2str(S1.w_grid(w_for_lambda))  ,').pdf' ]  );  close all;
end


%% Output Figures for Presentation
% 1. No bank-run
% clear all;  close all;
% load( [ '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model/Single State Solution/solution(par set 2).mat' ] );
% h = figure('position', [0, 0, 800, 600]);  linewidth=1;
% subplot(3,3,1);  plot( w_ode45, p_ode45,  'linewidth', linewidth );  title('price p', 'linewidth', linewidth);
% subplot(3,3,2); plot( w_ode45, sigmap_ode45 , 'linewidth', linewidth );  title('price volatility \sigma^p');
% subplot(3,3,3); plot( w_ode45, 0*psi_ode45, 'linewidth', linewidth );  title('price jump \kappa^p');
% subplot(3,3,4); plot( w_ode45, xg_ode45 , 'linewidth', linewidth);  title('liquidity premium x^g');
% subplot(3,3,5); plot( w_ode45, Rp_ode45 , 'linewidth', linewidth);  title('risk premium');
% subplot(3,3,6); plot( w_ode45(w_ode45>0.02), muw_ode45(w_ode45>0.02) , 'linewidth', linewidth);  title('growth rate of w');
% subplot(3,3,7); plot( w_ode45, SharpeRatio_ode45 , 'linewidth', linewidth);  title('Sharpe ratio');
% subplot(3,3,8); plot( w_ode45, SharpeRatio_ode45 , 'linewidth', linewidth);  title('adjusted Sharpe ratio');
% subplot(3,3,9); plot( w_ode45, r_ode45 , 'linewidth', linewidth);  title('interest rate');
% saveTightFigure( h, 'present_noJump.pdf' );
% 
% % 2. Constant intensity
% clear all;  close all;
% load( [ '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model/Single State Solution/solution(lambda=0.32).mat' ] );
% h = figure('position', [0, 0, 800, 600]);   linewidth=1;   
% subplot(3,3,1);  plot(  w_sol, p_new, w_ode45 , p_ode45 , ':m',   'linewidth', linewidth  );    title('price p');      legend( 'lambda=0.32' ,  'lambda=0'  );      
% subindex = w_sol>0.02 & solved_vec  ;   subindex_ode = w_ode45>0.02;
% subplot(3,3,2);  plot(  w_sol(subindex), sigmap_sol(subindex),  w_ode45(subindex_ode) , sigmap_ode45(subindex_ode) , ':m',   'linewidth', linewidth );    title('price volatility \sigma^p');     legend( 'lambda=0.32' ,  'lambda=0'  );      
% subindex = xK_sol<10&  w_sol~=0 & w_sol~=1 & solved_vec;  subindex_ode = xK_ode45<15& (psi_ode45<=1) & w_ode45~=0 & w_ode45~=1;
% subplot(3,3,3);  plot( w_sol,  kappap_sol,  w_sol ,  0*kappap_sol  , ':m' ,    'linewidth', linewidth );    title('price jump \kappa^p');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subplot(3,3,4);  plot( w_sol,  liq_prem_sol,  w_ode45,  zeros(N_ode45,1), ':m',   'linewidth', linewidth);    title('liquidity premium');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subplot(3,3,5);  plot( w_sol,  muR_sol + AH./p_sol,  w_ode45,   muR_ode45 + AH./p_ode45 , ':m',   'linewidth', linewidth);    title('risk premium');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subplot(3,3,6);  plot( w_sol,  muw_sol ,  w_ode45,   muw_ode45 , ':m',    'linewidth', linewidth);    title('growth rate of w');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subindex = w_sol>0.02 & solved_vec;
% subplot(3,3,7);  plot( w_sol(subindex), SharpeRatio(subindex) ,  w_ode45(subindex),   SharpeRatio_ode45(subindex) , ':m',  'linewidth', linewidth);    title('sharpe ratio');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subplot(3,3,8);  plot( w_sol(subindex),  AdjustedSR(subindex) ,  w_ode45(subindex),   SharpeRatio_ode45(subindex) , ':m',  'linewidth', linewidth);    title('adjusted sharpe ratio');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% subplot(3,3,9);  plot( w_sol(subindex),  r_sol(subindex) ,  w_ode45,   r_ode45 , ':m',  'linewidth', linewidth);    title('interest rate');    legend( 'lambda=0.32' ,  'lambda=0'  );      
% saveTightFigure( h,  [ 'present_ConstJump_lambda=', num2str(lambda), '.pdf' ]  );  

% 3. The whole model
clear all;  close all;
cd( '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model/Codes_2017-09-15/Results-Pivot' )
parset_choice = 2;
load( [ 'solution(par set ', num2str(parset_choice)  ,').mat' ] )
% Output w dimension
h = figure('position', [0, 0, 800, 600]);  linewidth=1;  
i1 = 10;  i2 = 5; 
subplot(3,3,1);  plot(  w_grid, p_2D_new(:,i1),   w_grid, p_2D_new(:,i2),  '--'   , w_ode45 , p_ode45 , ':m',   'linewidth', linewidth  );    title('price p');     xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],  '\lambda=0'  );
sub_w_index  =  w_grid>0.02 ;   subindex_ode = w_ode45>0.02;
subplot(3,3,2);  plot(  w_grid(sub_w_index), sigmap_matrix(sub_w_index, i1),  w_grid(sub_w_index), sigmap_matrix(sub_w_index, i2),  '--'  , w_ode45(subindex_ode) , sigmap_ode45(subindex_ode) , ':m',  'linewidth', linewidth );    title('price volatility \sigma^p');  xlabel('w');  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],  '\lambda=0'  );
subplot(3,3,3);  plot( w_grid(sub_w_index), kappap_matrix(sub_w_index, i1),  w_grid(sub_w_index), kappap_matrix(sub_w_index, i2),   '--'  ,  w_ode45 ,  0*N_ode45  , ':m' ,   'linewidth', linewidth );    title('price jump \kappa^p'); xlabel('w');  % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))],  '\lambda=0'  );
subplot(3,3,4);  plot( w_grid(sub_w_index), liq_prem_matrix(sub_w_index, i1),  w_grid(sub_w_index), liq_prem_matrix(sub_w_index, i2),  '--'  ,   w_ode45,  zeros(N_ode45,i1), ':m',  'linewidth', linewidth);    title('liquidity premium'); xlabel('w'); % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
subplot(3,3,5);  plot( w_grid(sub_w_index), RpH_matrix(sub_w_index, i1),  w_grid(sub_w_index), RpH_matrix(sub_w_index, i2),  '--',   w_ode45,   muR_ode45 + AH./p_ode45 , ':m',   'linewidth', linewidth);    title('risk premium'); xlabel('w'); % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
subindex_ode  =  muw_ode45<0.3 ;
subplot(3,3,6);  plot( w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i1),  w_grid(sub_w_index), muw_ratio_matrix(sub_w_index, i2),  '--'  , w_ode45(subindex_ode), muw_ratio_ode45(subindex_ode),  ':m',  'linewidth', linewidth);    title('growth rate of w'); xlabel('w');  % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
subplot(3,3,7);  plot( w_grid(sub_w_index), SharpeRatioH_matrix(sub_w_index, i1),  w_grid(sub_w_index), SharpeRatioH_matrix(sub_w_index, i2),  '--'  , w_ode45(sub_w_index), SharpeRatio_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('Sharpe ratio'); xlabel('w');  % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
subplot(3,3,8);  plot( w_grid(sub_w_index), SR_adjusted_matrix(sub_w_index, i1),  w_grid(sub_w_index), SR_adjusted_matrix(sub_w_index, i2),  '--'  ,  w_ode45(sub_w_index), SharpeRatio_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('adjusted Sharpe ratio'); xlabel('w');  % legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
subplot(3,3,9);  plot( w_grid(sub_w_index), r_matrix(sub_w_index, i1),  w_grid(sub_w_index), r_matrix(sub_w_index, i2),  '--'  ,  w_ode45(sub_w_index), r_ode45(sub_w_index),  ':m',  'linewidth', linewidth);    title('interest rate'); xlabel('w'); %  legend(   ['\lambda=', num2str(lambda_grid(i1))] ,   ['\lambda=', num2str(lambda_grid(i2))], '\lambda=0'  );
saveTightFigure( h,  [ 'present_w_dimension.pdf' ]  );  close all;
% Output lambda dimension
h = figure('position', [0, 0, 800, 600]);   linewidth=1;  
i1 = 4;      i2 = 15 ;  i3 = 25;
subplot(3,3,1);  plot(  [0;lambda_grid] , [ p_ode45(i1), p_2D_new(i1,:) ],   [0;lambda_grid],  [ p_ode45(i2), p_2D_new(i2,:)],  '--'  ,  [0;lambda_grid] , [ p_ode45(i3), p_2D_new(i3,:)] , ':',   'linewidth', linewidth  );    title('price p');     xlabel('lambda'); legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))],  ['w=', num2str(w_grid(i3))]  );
subplot(3,3,2);  plot(  [lambda_grid],  [sigmap_matrix(i1,:)]  ,    [lambda_grid],  [sigmap_matrix(i2,:)],  '--'  , [lambda_grid],  [sigmap_matrix(i3,:)], ':',  'linewidth', linewidth );    title('price volatility \sigma^p');  xlabel('lambda'); % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,3);  plot( [0;lambda_grid],  [ 0, kappap_matrix(i1,:) ] ,  [0;lambda_grid], [ 0,kappap_matrix(i2,:)],   '--'  ,[0;lambda_grid] ,  [0,kappap_matrix(i3,:)]  , ':' ,   'linewidth', linewidth );    title('price jump \kappa^p'); xlabel('lambda'); % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,4);  plot( [0;lambda_grid],   [0, liq_prem_matrix(i1,:)] ,  [0;lambda_grid],  [0, liq_prem_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [0,  liq_prem_matrix(i3,:)], ':',  'linewidth', linewidth);    title('liquidity premium'); xlabel('lambda');  % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,5);  plot( [0;lambda_grid],   [ Rp_ode45(i1), RpH_matrix(i1,:) ] ,  [0;lambda_grid], [ Rp_ode45(i2),RpH_matrix(i2,:)],  '--'  ,[0;lambda_grid],   [ Rp_ode45(i3),RpH_matrix(i3,:)] , ':',   'linewidth', linewidth);    title('risk premium'); xlabel('lambda');  % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,6);  plot( [ lambda_grid], [  muw_ratio_matrix(i1,:)] ,  [ lambda_grid],  [ muw_ratio_matrix(i2,:)],  '--'  ,[ lambda_grid],  [muw_ratio_matrix(i3,:)] ,':',  'linewidth', linewidth);    title('growth rate of w'); xlabel('lambda');  % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,7);  plot( [lambda_grid], [SharpeRatioH_matrix(i1,:)] ,  [lambda_grid], [SharpeRatioH_matrix(i2,:)],  '--'  ,[lambda_grid],  [SharpeRatioH_matrix(i3,:)] ,  'linewidth', linewidth);    title('Sharpe ratio'); xlabel('lambda');   % legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,8);  plot( [lambda_grid], [SR_adjusted_matrix(i1,:)] ,  [lambda_grid], [SR_adjusted_matrix(i2,:)],  '--'  ,[lambda_grid],  [SR_adjusted_matrix(i3,:)] ,':', 'linewidth', linewidth);    title('adjusted Sharpe ratio'); xlabel('lambda'); %  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
subplot(3,3,9);  plot( [0;lambda_grid], [r_ode45(i1), r_matrix(i1,:)] ,  [0;lambda_grid], [r_ode45(i2),r_matrix(i2,:)],  '--'  ,[0;lambda_grid],  [r_ode45(i3),r_matrix(i3,:)] , ':', 'linewidth', linewidth);    title('interest rate'); xlabel('lambda'); %  legend(   ['w=', num2str(w_grid(i1))] ,   ['w=', num2str(w_grid(i2))], ['w=', num2str(w_grid(i3))]  );
saveTightFigure( h,  [ 'present_lambda_dimension.pdf' ]  );  close all;

% 4. Comparison of different Q
load('solution(par set 1).mat'); 
S1=load('solution(par set 1).mat');     S2=load('solution(par set 3).mat');  
labelS1_1 = ['Q=', num2str(S1.multiplier)];  labelS1_2= ['Q=', num2str(S1.multiplier), ' no jump'];    
labelS2_1 = ['Q=', num2str(S2.multiplier)];  labelS2_2= ['Q=', num2str(S2.multiplier), ' no jump'];    
if( option_Q == 4 )
    labelS1_1 = ['Countercyclical Q'];  labelS1_2= ['Countercyclical Q, no jump'];    
    labelS2_1 = ['Cyclical Q'];  labelS2_2= ['Cyclical Q, no jump'];  
end

%% w dimension
for( index_lambda = [4 11]  )
index_w  =  w_grid>0.02 ;   index_ode=w_ode45>0.02;
h = figure('position', [0, 0, 800, 600]);  linewidth=1;  
% Part 1: price comparison
subplot(3,3,1);
plot(  S1.w_grid , S1.p_2D_new(:,index_lambda),  '--',  'linewidth', linewidth  );    title('price p');     xlabel('w');  hold on;
plot(  S2.w_grid , S2.p_2D_new(:,index_lambda),  'linewidth', linewidth  );    title('price p');     xlabel('w'); 

% Part 2: sigmap comparison
subplot(3,3,2);
plot(  S1.w_grid(index_w) , S1.sigmap_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('price volatility \sigma^p');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.sigmap_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('price volatility \sigma^p');     xlabel('w'); 
legend( labelS1_1, labelS2_1 ); 

% Part 3: kappap comparison
subplot(3,3,3);
plot(  S1.w_grid(index_w) , S1.kappap_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('price jump \kappa^p');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.kappap_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('price jump \kappa^p');     xlabel('w'); 
   
% Part 4: liquidity premium comparison
subplot(3,3,4);
plot(  S1.w_grid(index_w) , S1.liq_prem_matrix(index_w,index_lambda) , '--',  'linewidth', linewidth  );     title('liquidity premium');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.liq_prem_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('liquidity premium');     xlabel('w'); 
   
% Part 5: risk premium comparison
subplot(3,3,5);
plot(  S1.w_grid(index_w) , S1.RpH_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('risk premium');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.RpH_matrix(index_w,index_lambda), 'linewidth', linewidth  );   title('risk premium');     xlabel('w'); 
   
% Part 6 : growth rate comparison
subplot(3,3,6);  index_ode = w_ode45>0.01;
plot(  S1.w_grid(index_w) , S1.muw_ratio_matrix(index_w,index_lambda), '--',  'linewidth', linewidth  );     title('growth rate of w');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.muw_ratio_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('growth rate of w');     xlabel('w'); 
   
% Part 7: Sharpe Ratio comparison
subplot(3,3,7);
plot(  S1.w_grid(index_w) , S1.SharpeRatioH_matrix(index_w,index_lambda), '--', 'linewidth', linewidth  );     title('Sharpe Ratio');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.SharpeRatioH_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('Sharpe Ratio');     xlabel('w'); 
   
% Part 8: adjusted Sharpe Ratio comparison
subplot(3,3,8);
plot(  S1.w_grid(index_w) ,  max(S1.SR_adjusted_matrix(index_w,index_lambda), 0) , '--',  'linewidth', linewidth  );     title('adjusted Sharpe Ratio');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) ,  max(S2.SR_adjusted_matrix(index_w,index_lambda),0) ,  'linewidth', linewidth  );   title('adjusted Sharpe Ratio');     xlabel('w'); 
   
% Part 9: interet rate comparison
subplot(3,3,9);
plot(  S1.w_grid(index_w) , S1.r_matrix(index_w,index_lambda), '--', 'linewidth', linewidth  );     title('interest rate');  xlabel('w');  hold on;
plot(  S2.w_grid(index_w) , S2.r_matrix(index_w,index_lambda),  'linewidth', linewidth  );   title('interest rate');     xlabel('w'); 

saveTightFigure( h,  [ 'present_LiquidityEffects_w(lambda=',  num2str(  S1.lambda_grid( index_lambda )  )  ,').pdf' ]  );  close all;

end

%% lambda dimension
for( w_for_lambda = [6 27 40]  )
lambda_index = lambda_grid>=0.1;  
h = figure('position', [0, 0, 800, 600]);  linewidth=1;  
% Part 1: price comparison
subplot(3,3,1);
plot(  [0;S1.lambda_grid] , [ S1.p_ode45(w_for_lambda), S1.p_2D_new(w_for_lambda,:)]  ,  '--', 'linewidth', linewidth  );    title('price p');     xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [ S2.p_ode45(w_for_lambda), S2.p_2D_new(w_for_lambda,:)],  'linewidth', linewidth  );    title('price p');     xlabel('\lambda'); 
 
% Part 2: sigmap comparison
subplot(3,3,2);
plot(  [0;S1.lambda_grid(lambda_index)  ] ,  [S1.sigmap_ode45(w_for_lambda),S1.sigmap_matrix(w_for_lambda,lambda_index)] , '--', 'linewidth', linewidth  );     title('price volatility \sigma^p');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid(lambda_index) ] ,  [S2.sigmap_ode45(w_for_lambda),S2.sigmap_matrix(w_for_lambda,lambda_index)],  'linewidth', linewidth  );   title('price volatility \sigma^p');     xlabel('\lambda'); 
legend( labelS1_1,  labelS2_1 );
 
% Part 3: sigmap comparison
subplot(3,3,3);
plot(  [0;S1.lambda_grid] , [0,S1.kappap_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('price jump \kappa^p');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [0,S2.kappap_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('price jump \kappa^p');     xlabel('\lambda'); 
 
% Part 4: liquidity premium comparison
subplot(3,3,4);
plot(  [0;S1.lambda_grid] , [0,S1.liq_prem_matrix(w_for_lambda,:)] , '--',  'linewidth', linewidth  );     title('liquidity premium');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [0,S2.liq_prem_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('liquidity premium');     xlabel('\lambda'); 
 
% Part 5: risk premium comparison
subplot(3,3,5);
plot(  [0;S1.lambda_grid] , [S1.Rp_ode45(w_for_lambda),S1.RpH_matrix(w_for_lambda,:)], '--', 'linewidth', linewidth  );     title('risk premium');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.Rp_ode45(w_for_lambda),S2.RpH_matrix(w_for_lambda,:)],  'linewidth', linewidth  );   title('risk premium');     xlabel('\lambda'); 
 
% Part 6: growth rate of w
subplot(3,3,6);
plot(  [S1.lambda_grid] , [S1.muw_ratio_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('growth rate of w');  xlabel('\lambda');  hold on;
plot(  [S2.lambda_grid] , [S2.muw_ratio_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('growth rate of w');     xlabel('\lambda'); 
 
% Part 7: Sharpe Ratio comparison
subplot(3,3,7);
plot(  [0;S1.lambda_grid] , [S1.SharpeRatio_ode45(w_for_lambda),S1.SharpeRatioH_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('Sharpe Ratio');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.SharpeRatio_ode45(w_for_lambda),S2.SharpeRatioH_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('Sharpe Ratio');     xlabel('\lambda'); 
 
% Part 8: adjusted Sharpe Ratio comparison
subplot(3,3,8);
plot(  [0;S1.lambda_grid] , [S1.SharpeRatio_ode45(w_for_lambda),S1.SR_adjusted_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('adjusted Sharpe Ratio');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.SharpeRatio_ode45(w_for_lambda),S2.SR_adjusted_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('adjusted Sharpe Ratio');     xlabel('\lambda'); 
 
% Part 9: interest rate
subplot(3,3,9);
plot(  [0;S1.lambda_grid] , [S1.r_ode45(w_for_lambda),S1.r_matrix(w_for_lambda,:)], '--',  'linewidth', linewidth  );     title('interest rate');  xlabel('\lambda');  hold on;
plot(  [0;S2.lambda_grid] , [S2.r_ode45(w_for_lambda),S2.r_matrix(w_for_lambda,:)], 'linewidth', linewidth  );   title('interest rate');     xlabel('\lambda'); 
 
saveTightFigure( h,  [ 'present_LiquidityEffects_lambda(w=',  num2str(  S1.w_grid( w_for_lambda )  )  ,').pdf' ]  );  close all;
end


end


%% Cross Q-Regime comparison. 
MainPath = '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model/Codes_2017-09-15';

cd( MainPath )
SolutionName = 'solution(par set 3).mat';
S1= load([MainPath,'/Results-Constant Q/',SolutionName]);   
S2 =load([MainPath,'/Results-Linear Q/',SolutionName]);   
S3 =  load([MainPath,'/Results-Quadratic Q/',SolutionName]);   

figure_size = [0, 0, 600, 400]; linewidth=2;  
figure_size2 = [0, 0, 800, 400];
w_subgrid = w_grid(w_grid<0.6);   w_grid_index = find(w_grid<0.6);
lambda_ind = 8;    w_ind = 15;
legend_text = { ['Q=', num2str(S1.multiplier)] ,   ['Q=', num2str(1.4*S2.multiplier),'(1-w)'], ['Q=', num2str(2*S3.multiplier),'(1-w)^2'] } ;
xlab1 = ['w dimension at lambda=', num2str(lambda_grid(lambda_ind))];  % For w dimension
xlab2 = ['\lambda dimension at w=', num2str(w_grid(w_ind))] ;  % For lambda dimension

% A comparison of total supply of government bonds
h = figure('position', figure_size);  
plot( w_subgrid  ,   S1.Q(w_subgrid) ,    w_subgrid,    S2.Q(w_subgrid) ,  '--',  w_subgrid,  S3.Q(w_subgrid), '-.' )
title('Liquidity Supply'); xlabel('w');  legend(  legend_text );
saveTightFigure( h,  [ 'Q_Comparison_Qfunction.pdf' ]  );  close all;

% Compare the price effects
h = figure('position', figure_size2); 
subplot(1,2,1); plot( w_grid  ,   S1.p_matrix(:,lambda_ind) ,    w_grid,    S2.p_matrix(:,lambda_ind)  ,  '--',  w_grid,  S3.p_matrix(:,lambda_ind) , '-.' )
title('Price'); xlabel( xlab1 );  legend( legend_text  );
subplot(1,2,2); plot( lambda_grid  ,   S1.p_matrix(w_ind,:) ,    lambda_grid,    S2.p_matrix(w_ind,:) ,  '--',  lambda_grid,  S3.p_matrix(w_ind,:), '-.' )
title('Price'); xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_Price.pdf' ]  );  close all;

% Compare the recovery speed
h = figure('position', figure_size2); 
SLM1 = slmengine( w_grid ,  S1.muw_ratio_matrix(:,lambda_ind)  ,'plot', 'off',  'concavedown', 'on', 'knots',   floor(numel(w_grid)/4)   );  
SLM2 =  slmengine( w_grid ,  S2.muw_ratio_matrix(:,lambda_ind)  ,'plot', 'off',  'concavedown', 'on', 'knots',   floor(numel(w_grid)/4)   );  
SLM3 =  slmengine( w_grid ,  S3.muw_ratio_matrix(:,lambda_ind)  ,'plot', 'off',  'concavedown', 'on', 'knots',   floor(numel(w_grid)/4)   );  
subplot(1,2,1);  plot( w_grid  ,   slmeval( w_grid, SLM1 ) ,    w_grid,     slmeval( w_grid, SLM2 ) ,  '--',  w_grid,   slmeval( w_grid, SLM3 ) , '-.' )
title('Growth rate \mu^w'); xlabel( xlab1 );  legend( legend_text  );
w_ind = 5;
subplot(1,2,2);  plot( lambda_grid  ,   S1.muw_ratio_matrix(w_ind,:) ,    lambda_grid,    S2.muw_ratio_matrix(w_ind,:)  ,  '--',  lambda_grid,  S3.muw_ratio_matrix(w_ind,:) , '-.' )
title('Growth rate \mu^w'); xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_GrowthRatew.pdf' ]  );  close all;

% Compare the interest rate
h = figure('position', figure_size2);   
w_ind = 15;
subplot(1,2,1);  plot( w_grid  ,   S1.r_matrix(:,lambda_ind) ,    w_grid,    S2.r_matrix(:,lambda_ind)  ,  '--',  w_grid,  S3.r_matrix(:,lambda_ind) , '-.' )
title('Interest rate'); xlabel( xlab1 );  legend( legend_text  );
subplot(1,2,2);  plot( lambda_grid  ,   S1.r_matrix(w_ind,:) ,    lambda_grid,    S2.r_matrix(w_ind,:)  ,  '--',  lambda_grid,  S3.r_matrix(w_ind,:) , '-.' )
title('Interest rate'); xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_InterestRate.pdf' ]  );  close all;

% Compare liquidity premium
h = figure('position', figure_size2);   
subplot(1,2,1);  plot( w_grid  ,   S1.liq_prem_matrix(:,5) ,    w_grid,    S2.liq_prem_matrix(:,5)  ,  '--',  w_grid,  S3.liq_prem_matrix(:,5) , '-.' )
title('Liquidity Premium');  xlabel( ['w dimension at lambda=', num2str(lambda_grid(5))] );  legend( legend_text  );
subplot(1,2,2);  plot( [0; lambda_grid]  ,   [0, S1.liq_prem_matrix(w_ind,:)] ,     [0;lambda_grid],    [0,S2.liq_prem_matrix(w_ind,:)]  ,  '--',  [0; lambda_grid],  [0,S3.liq_prem_matrix(w_ind,:)] , '-.' )
title('Liquidity Premium');  xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_LiqPrem.pdf' ]  );  close all;

% Compare risk premium
h = figure('position', figure_size);   
subplot(1,2,1);  plot( w_grid  ,   S1.RpH_matrix(:,lambda_ind)  ,    w_grid,    S2.RpH_matrix(:,lambda_ind)   ,  '--',  w_grid,  S3.RpH_matrix(:,lambda_ind) , '-.' )
title('Risk Premium');  xlabel( xlab1 );
subplot(1,2,2);  plot( lambda_grid  ,   S1.RpH_matrix(w_ind,:) ,    lambda_grid,    S2.RpH_matrix(w_ind,:)  ,  '--',  lambda_grid,  S3.RpH_matrix(w_ind,:) , '-.' )
title('Risk Premium');  xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_RiskPrem.pdf' ]  );  close all;

% Compare the disruption
h = figure('position', figure_size2);   
lambda_ind = 8;
subplot(1,2,1); plot( w_grid  ,   S1.kappap_matrix(:,lambda_ind) ,    w_grid,    S2.kappap_matrix(:,lambda_ind)  ,  '--',  w_grid,  S3.kappap_matrix(:,lambda_ind) , '-.' )
title('Price Jump \kappa^p'); xlabel( xlab1 );  legend( legend_text  );
w_ind = 15;
subplot(1,2,2); plot( lambda_grid  ,   S1.kappap_matrix(w_ind,:) ,    lambda_grid,    S2.kappap_matrix(w_ind,:)  ,  '--',  lambda_grid,  S3.kappap_matrix(w_ind,:) , '-.' )
title('Price Jump \kappa^p'); xlabel( xlab2 );  legend( legend_text  );
saveTightFigure( h,  [ 'Q_Comparison_Jump.pdf' ]  );  close all;





